clip-02-05

clip-02-05.jl

Load Julia packages (libraries) needed for the snippets in chapter 0

using StatisticalRethinking, Optim
gr(size=(600,300));
Plots.GRBackend()

snippet 3.2

Grid of 1001 steps

p_grid = range(0, step=0.001, stop=1);
0.0:0.001:1.0

all priors = 1.0

prior = ones(length(p_grid));
1001-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 ⋮
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

Binomial pdf

likelihood = [pdf(Binomial(9, p), 6) for p in p_grid];
1001-element Array{Float64,1}:
 0.0
 8.374825191600018e-17
 5.3438084689919926e-15
 6.068652771862778e-14
 3.399517250519025e-13
 1.2929107734374954e-12
 3.848982544705552e-12
 9.676432504149003e-12
 2.1495830280142858e-11
 4.344655104237097e-11
 ⋮
 4.098446591204723e-5
 2.7622876204442173e-5
 1.7500535729793716e-5
 1.0188911348240822e-5
 5.248259379330843e-6
 2.227480958032322e-6
 6.639762126411521e-7
 8.34972583212597e-8
 0.0

As Uniform prior has been used, unstandardized posterior is equal to likelihood

posterior = likelihood .* prior;
1001-element Array{Float64,1}:
 0.0
 8.374825191600018e-17
 5.3438084689919926e-15
 6.068652771862778e-14
 3.399517250519025e-13
 1.2929107734374954e-12
 3.848982544705552e-12
 9.676432504149003e-12
 2.1495830280142858e-11
 4.344655104237097e-11
 ⋮
 4.098446591204723e-5
 2.7622876204442173e-5
 1.7500535729793716e-5
 1.0188911348240822e-5
 5.248259379330843e-6
 2.227480958032322e-6
 6.639762126411521e-7
 8.34972583212597e-8
 0.0

Scale posterior such that they become probabilities

posterior = posterior / sum(posterior);
1001-element Array{Float64,1}:
 0.0
 8.374825191541396e-19
 5.3438084689545866e-17
 6.068652771820298e-16
 3.3995172504952293e-15
 1.2929107734284453e-14
 3.84898254467861e-14
 9.67643250408127e-14
 2.149583027999239e-13
 4.3446551042066853e-13
 ⋮
 4.0984465911760347e-7
 2.7622876204248817e-7
 1.7500535729671214e-7
 1.01889113481695e-7
 5.248259379294106e-8
 2.22748095801673e-8
 6.639762126365043e-9
 8.349725832067523e-10
 0.0

snippet 3.3

Sample using the computed posterior values as weights

N = 10000
samples = sample(p_grid, Weights(posterior), N);
10000-element Array{Float64,1}:
 0.863
 0.387
 0.663
 0.443
 0.697
 0.646
 0.71
 0.818
 0.782
 0.762
 ⋮
 0.898
 0.609
 0.704
 0.42
 0.456
 0.743
 0.738
 0.723
 0.672

In StatisticalRethinkingJulia samples will always be stored in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), [:toss]);
Object of type Chains, with data of type 10000×1×1 Array{Union{Missing, Float64},3}

Log evidence      = 0.0
Iterations        = 1:10000
Thinning interval = 1
Chains            = Chain1
Samples per chain = 10000
parameters        = toss

parameters
      Mean    SD   Naive SE  MCSE    ESS
toss 0.6352 0.1387   0.0014 0.0014 9570.62

Describe the chain

describe(chn)
Log evidence      = 0.0
Iterations        = 1:10000
Thinning interval = 1
Chains            = Chain1
Samples per chain = 10000
parameters        = toss

Empirical Posterior Estimates:
==========================================
parameters
      Mean    SD   Naive SE  MCSE    ESS
toss 0.6352 0.1387   0.0014 0.0014 9570.62

Quantiles:
==========================================
parameters
      2.5% 25.0% 50.0% 75.0% 97.5%
toss 0.098 0.542 0.643 0.739 0.971

Plot the chain

plot(chn)
0 2500 5000 7500 10000 0.2 0.4 0.6 0.8 toss Iteration Sample value 0.00 0.25 0.50 0.75 1.00 0.0 0.5 1.0 1.5 2.0 2.5 toss Sample value Density

Compute the MAP (maximumaposteriori) estimate

x0 = [0.5]
lower = [0.0]
upper = [1.0]

function loglik(x)
  ll = 0.0
  ll += log.(pdf.(Beta(1, 1), x[1]))
  ll += sum(log.(pdf.(Binomial(9, x[1]), repeat([6], N))))
  -ll
end

(qmap, opt) = quap(samples, loglik, lower, upper, x0)
([0.666667], Results of Optimization Algorithm
 * Algorithm: Fminbox with Gradient Descent
 * Starting Point: [0.5]
 * Minimizer: [0.6666666666544907]
 * Minimum: 1.297811e+04
 * Iterations: 4
 * Convergence: true
   * |x - x'| ≤ 0.0e+00: false
     |x - x'| = 5.74e-11
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: true
     |f(x) - f(x')| = 0.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: false
     |g(x)| = 1.80e-06
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 108
 * Gradient Calls: 108)

Show optimization results

opt
Results of Optimization Algorithm
 * Algorithm: Fminbox with Gradient Descent
 * Starting Point: [0.5]
 * Minimizer: [0.6666666666544907]
 * Minimum: 1.297811e+04
 * Iterations: 4
 * Convergence: true
   * |x - x'| ≤ 0.0e+00: false
     |x - x'| = 5.74e-11
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: true
     |f(x) - f(x')| = 0.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: false
     |g(x)| = 1.80e-06
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 108
 * Gradient Calls: 108

Fit quadratic approcimation

quapfit = [qmap[1], std(samples, mean=qmap[1])]
2-element Array{Float64,1}:
 0.6666666666544907
 0.14222160504709933

snippet 3.4

Create a vector to hold the plots so we can later combine them

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")
0 2500 5000 7500 10000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Draws

snippet 3.5

Analytical calculation

w = 6
n = 9
x = 0:0.01:1
p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
0.00 0.25 0.50 0.75 1.00 0 1 2 3 4 5 Sample density Conjugate solution

Add quadratic approximation

plot!( p[2], x, pdf.(Normal( quapfit[1], quapfit[2] ) , x ), lab="Quap approximation")
plot(p..., layout=(1, 2))
0 2500 5000 7500 10000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Draws 0.00 0.25 0.50 0.75 1.00 0 1 2 3 4 5 Sample density Conjugate solution Quap approximation

End of clip_02_05.jl

This page was generated using Literate.jl.